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Abstract 

We study efficiency of higher order integrator schemes for the hybrid Monte Carlo 
(HMC) algorithm. Numerical tests are performed for Quantum Chromo Dynamics 
(QCD) with two flavors of Wilson fermions. We compare 2nd, 4th and 6th order 
integrators at various quark masses. The performance depends on both volume 
and quark mass. On currently accessible large lattices ( V ~ 24 4 ), higher order 
integrators can be more efficient than the 2nd order one only in heavy quark region, 
m q a > 0.3. Thus we conclude that for most full QCD simulations, except for heavy 
quark case, the usual 2nd order integrator is the best choice. 



1 Introduction 



Inclusion of dynamical fermions is one of major difficulties in lattice QCD simula- 
tions since eventually one finds that the simulations require huge computational time. 
The standard algorithm for full QCD simulations is the hybrid Monte Carlo (HMC) 
algorithm 0]. While the basic idea of the HMC is a combination of molecular dynamics ( 
MD ) and Metropolis accept/reject steps, performance of its algorithm depends on tun- 
ing: matrix solver, parameter tuning, integration scheme etc. The matrix solver appears 
in the fermionic force calculations and takes the dominant time in the HMC simulations. 
The choice of an efficient matrix solver is an important subject to reduce CPU time [Q]. 
Parameters ( (3 and k etc ) of a MD Hamiltonian can be tuned so that the Metropolis 
acceptance rate increases ||]. 

One may choose any integrator for the MD step provided that the following two 
conditions are satisfied: 

• area preserving 

• time reversibility 

Usually the ( 2nd order ) leapfrog integrator is used for the HMC. The leading 
integration errors are 0(At 3 ), where At is the step size of an elementary MD step. Due 
to these errors the Hamiltonian is not conserved. Let AH be an energy violation at the 
end of a MD trajectory. To achieve the correct equilibrium a new configuration should 
be accepted by a global Metropolis test with a probability: 

P oc min(l, exp(-AH)). (1) 

In order to have a high acceptance we may consider a more accurate integration 
scheme to reduce AH. The multiple time scale method @] which removes the dominant 
errors from the gauge part worked well. An idea which controls the integration errors 
with an adaptive step size was also explored. However no practical gain appeared for 
QCD case Q. One may employ a higher order integrator which has higher order integra- 
tion errors in At. In general higher order integrators need more arithmetic operations 
than the 2nd order one. Therefore it is non-trivial whether the higher order integrators 
serve as an efficient speed-up source to the HMC. For QCD the 4th order one was studied 
on a 4 4 lattice §, and was found not to be efficient enough on such a small lattice. 

When one compares higher order integrators, volume dependence must be considered. 
The average acceptance of an n-th order integrator is given by ~ erfc(cU 1 / 2 At n )|, where 
V is volume of the system considered and c is a constant. To keep a constant acceptance, 
At should scale ~ l/U 1//2n . This scaling behaviour suggests that the higher order one 
will be efficient for a lattice bigger than a certain size. However we do not know the 
value of this lattice size. In this study we perform HMC simulations with 2nd, 4th and 
6th order integrators and clarify which integrator is efficient for a given lattice. It now 
becomes feasible to perform a simulation on rather big lattices as 24 3 x 40 — 48 lattices^. 
So it is worthwhile to study whether, on such lattices, the higher order integrators are 
more efficient than the standard leapfrog ( 2nd ) one. 

'See Sec. 4. 
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In Sec. 2 we describe the lattice QCD action used in our HMC simulations. In Sec. 3 
we describe the higher order integrators which we use. In Sec. 4 we discuss the optimal 
efficiency, acceptance and step size of the HMC. In Sec. 5 we give a criterion to compare 
various integrators. In Sec. 6 we present our numerical results. Finally we summarize our 
results in Sec. 7. 

2 Lattice QCD action 

We use the standard plaquette gauge action and two flavors Wilson fermion action ||. 
The partition function is given by 

Z = JvUdet[M(U)*M(U)]exp(-Sg), (2) 

where U stands for SU(3) link variables and 

S 9 = ^Tr(l-U p ) (3) 

where U p stands for the plaquette and j3 is the gauge coupling, and the Wilson fermion 
matrix M(U) is given by 

M^U) = Sij + ^1(1, ~ l)tfi,Ai+^ - (7m + l)Ul^S i>j+l ,] (4) 

where n is the hopping parameter. 

The expectation value of some operator 0(U) is given by 

< o >= J VUO(U) det[M(U^M(U)} exp(-S g )/Z. (5) 

Using pseudofermion fields 4> we replace the determinant in Eq.(j2|) with a path- 
integral as 

Z = j VUVfV(j)eK-p[-Sg- <f) ] {M(U) ] M(U)y l <i)]. (6) 

Introducing momenta p conjugate to link variables we define the Hamiltonian used 
in the HMC as 

H = \p 2 + S g + ^{M(U)Hl(U))- V (7) 
and the partition function will be 

Z = J VUV0*V0Vpexp(-H), (8) 

which gives the same expectation values as that from Eq.(Q). 
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3 Higher order integrators 



In this section we describe higher order integrator scheme which we use for the present 
study. Let H be a classical Hamiltonian, 

H = \ P 2 + S(q) (9) 

where q = (q\, q2, ■ ■■) and p = (pi,P2> •••) are coordinate variables and conjugate momenta 
respectively, and S(q) represents a potential term of the system. For simplicity we use 
scaler variables p and q. The same discussion applies for QCD case where SU(3) link 
variables are used. 

In the MD step we solve Hamilton's equations, 

dq i = dH 

dt d Pi y ' 

dpi dH , . 

— = fill 

dt d qi ' { ' 

approximately by an appropriate integrator. In general these equations are not solvable 
analytically. Let T^p^At) be an elementary MD step with a time interval ( step size ) 
At, which evolves (p, q) to (p',q')- 

T MD {At) : (p,q) — {p',q'). (12) 

Requirements of the HMC to the integrator Tj\^£)(Ai) are (a) time reversible: 

T MD (-At):(p',q')-^{p,q) (13) 

and (b) area preserving: 

dpdq = dp'dq' (14) 

i.e. invariance of the measure. All integrators having the above requirements can be 
used for the HMC. The simplest integrator is the 2nd order leapfrog method which has 
been commonly used in the HMC of the current full QCD simulations. The 2nd order 
leapfrog scheme is explicitly written as 



q(t + f) =q(t) + fp(t) 

p{t + At) =p(t)-At d J^^l (15) 
q(t + At) =q(t + f) + fp(t + At). 

While we start the integrator with variables q, alternatively we can use momenta p for the 
starting variables. This 2nd order leapfrog integrator causes C(At 3 ) integration error. 
In order to construct a class of higher order integrators, it is convenient to use the 



Lie algebraic formalism J|, 11, |12|, 14]. The Hamilton's equation is written as 



! = {/,*> as) 
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where / = p or q, and {,} stands for the Poisson bracket, i.e. 

{/ ,* } = B v«E (17) 

Y ^ 5 K d H 9 Pi 

Defining the linear ( Lie ) operator L(H) as 

L(H)f = {f,H}, (18) 
we have the formal solution of the Hamilton's equation: 

f(t + At) = exp(AtL(H))f(t). (19) 
Since L(-) is a linear operator, we have 

L{H) = L(\p 2 ) + L(S(q)) (20) 
= T + V (21) 

where T = L(^p 2 ) and V = L(S(q)) stand for a kinetic and potential terms respectively. 
Using the Lie algebraic formalism, one finds that the 2nd leapfrog integrator corresponds 
to a decomposition of the exponential in Eq,(|i"9D as 

/(t + At) = exp(At(T + V))/(t) 

= {exp(^AtT)exp(AtV)exp(-AtT) + 0(At 3 )}f(t). (22) 

Note that T and V do not commute with each other and 0(At 3 ) decomposition errors 
appear. An important observation here is that the order of the decomposition error 
coincides with that of the integration error. 

An arbitrary order integrator can be found by decomposing the exponential with the 
desired order. In general, the exponential is decomposed as 

exp(At(T + V))=Y[ exp{ciAtT)exp{diAtV) + 0(At n+1 ) (23) 

i 

where q and di are determined so that the decomposition is correct up to 0(At n ). It 
is not obvious how to obtain such Cj and di in any order. Fortunately higher ei>en-order 



integrators are known to be constructed from a combination of lower order integrators! 11 



12, 13]. Let us call Cr2nd(Ai) the 2nd order decomposition ( or integrator), 



G 2nd (Ai) = exp(^AtT)exp(AtV)exp(^AtT). (24) 



The 4th order integrator is given by a product of three 2nd order integrators [10, 11, 
G 4th (At) = G 2 nd{aiAt)G2nd{a2^t)G 2 nd{aiAt) (25) 



where the coefficients aj are given by 



ai = 232178' ^ 
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2 l/3 

'2-2V3- 



(27) 



This construction scheme is easily generalized for an arbitrary even-order one[ll, [l2|, [II 
(2k+2)-th order integrator is given recursively by 



G 2k+2 (At) = G 2fc (&iAt)G 2fe (6 2 At)G 2fc (6iAt), 
where the coefficient bi are 

h 



1 



2 _ 2V(2*+i) 

2 1/(2A+1) 



(28) 



(29) 



(30) 



z 2 - 2 1 /(2fc+l) ' 

While one can find an arbitrary higher even-order integrator with Eq. , the number 
of elementary steps ( 2nd order integrator ) grows with the order of the integrator as 
follows: 

f 1 2nd 



# of 2nd order integrator 
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4th 
6th 



(31) 



3 n / 2 ~ l nth. 



A bottleneck of the HMC for QCD is the force calculation dS/dq which needs a large 
amount of computational time devoted to a matrix solver. Except for some small over- 
head, the computational cost of the HMC is proportional to the number of the force 
calculations. The 2nd order integrator contains one force calculation. Therefore the 
computational cost of the higher order algorithm can be counted by Eq.(^Tj), which 
indicates that the cost grows rapidly with the order. 

If one can find a higher order integrator consisting of fewer force calculations it may 



be useful for HMC. In Ref[12|, such a higher order scheme is found numerically. In this 
study we also use the 6th order integrators of Ref|I^] which consist of seven 2nd order 
integrators instead of nine as in Eq. (|3l|) . The 6th order integrators of Ref|12] are written 
as 



G uh (At) = G 2nd (w 3 At)G 2nd (w 2 At)G 2nd (w 1 At)G 2nd (w At) 
x G 2nd (wi At ) G 2nd (w 2 At)G 2nd (w 3 At) , 



(32) 



where values of (wq, W\, w 2 , W3) are listed in Table 1. 

Note that all the higher even-order integrators described here satisfy the time re- 
versible and area preserving conditions since those are a product of the 2nd order inte- 
grators having the area preserving condition and are constructed in a symmetric way ( 
G(At)G(—At) = 1 ) which yield the time reversible condition. 



4 Optimal efficiency, acceptance and step size 

In this section we introduce an efficiency function which characterizes the speed of al- 
gorithm and derive formulae for the optimal acceptance and step size which define the 
optimal efficiency. 
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Yl 


Y2 


Y3 




-0.117767998417887e-l 


-0.2132285222000144e+l 


0.152886228424922e-2 


W 2 


0.235573213359357e+0 


0.426068187079180e-2 


-0.214403531630539e+l 


W 3 


0.784513610477560e+0 


0.143984816797678e+l 


0.144778256239930e+l 



Table 1: Parameter sets (Y1-Y3) of the 6th order integrators by Yoshida [12]. u>o is given 
by wq = 1 — w\ — u>2 — W3 . 



We define the efficiency function Eff by a product of step size At and acceptance 

Pacc(At): 

Eff (At) = P acc (At)At. (33) 

High Eff results in fast Markov step when producing configurations. A particular length 
of trajectory does not affect Eff since the acceptance stays almost constant for any 



trajectory length longer than a certain characteristic length [15]. This is a feature of the 



symplectic type integrator [14]. We fix the trajectory length to the unit length (=1) 
which is sufficiently longer than the characteristic length for the present studyQ. 

The acceptance decreases as At increases. In both limits of At = and oo, Eff goes 
to zero. We expect that Eff has one maximum at a certain At, which we call optimal 
step size At opt . 

Using At opt we define the optimal acceptance (P a cc)opt- 

(Pace) opt — Pacc(At Q pt) , (34) 

and the optimal efficiency (Eff) opt : 

(Eff) opt = m&x(Eff(At))= Eff (At opt ) (35) 

= (Pace) opt Atopt- (36) 

When we have (Eff) op t the maximum speed of the algorithm is achieved. Comparison 
among various integrators will be done with this (Eff) op t- At this stage, however, it is 
not obvious how to obtain (Eff) op t easily from Monte Carlo (MC) simulations. In the 
following consideration, we show that one coefficient governs (Eff) op t and it is easily 
obtainable from a MC simulation. 

At large volumes, the average acceptance with an average energy difference is evalu- 
ated as jn| 

(P acc ) = erfc(i(A^) 1 /2). (37) 

Although Eq. ([}?]) is applicable for the whole range of the acceptance, i.e. 1 > (P aC c) > 0, 
we are not interested in low acceptance with which the algorithm may not be efficient. 
Typically we may need (P a cc) > 50%. Instead of Eq.fl37|) we propose a simple exponential 
type formula: 

(P acc ) = eM-^(^H 2 ) 1/2 )- (38) 

^Several simulations have been done with both trajectory lengths of 0.5 and 1.0. We found no 
significant change in the acceptance among them. 
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In the limit of (AH) — > 0, Eq.(|38[) coincides with Eq.(p7D ( See Appendix ). Even if 
(AH) is not small enough ( (AH) ~ 1 ), from numerical tests we confirm that Eq,(|38|) 
is a good approximation to the exact value. Fig.l shows comparison of MC results 
and values from Eq, (|38|) . Fig. 2 shows a normalized error [(MC— Eq.(j38|)) /Eq.([38|)l as a 
function of (AH 2 ) 1 / 2 . From the comparison of MC results and Eq.(|38[) we notice that 
Eq.(|3q) agrees quite well with the MC results within 5% error up to (AH 2 ) 1 / 2 < 3, which 
roughly corresponds to (P aC c) > 20%. Therefore we adopt Eq.(|38|) as our formula for 

(Pace) • 

Let us now discuss At dependence of (AH 2 ) 1 / 2 which appears in Eq. (|38|) . The 2nd 
order integrator causes 0(At 3 ) integration errors after an elementary MD step. From 
the discussion of Ref . |l5|l , however, AH ~ At 2 rather than ~ At 3 . This comes from the 
fact that AH does not grow with the trajectory length ( = the number of elementary 
MD steps xAt ). Using (AH) ^(AH 2 ) which is expected from Creutz's equality fill 
(exp(-AH)) = 1 at small AH, we obtain 

(AH) ~ (AH 2 ) ~ VAt A (39) 

where V is the volume of the system. When we apply the same discussion for the n-th 
order integrator we have 

(AH 2 ) ~ VAt 2n . (40) 

Thus, 

(AH 2 ) 1 / 2 ~ V 1/2 At n . (41) 
For our later use we rewrite Eq.(^Tj) as 

(AH 2 ) 1 / 2 = C n V x / 2 At n + 0(At n+1 ), (42) 

where C n is a Hamiltonian ( model ) dependent coefficient, which is not known a priori. 



We checked Eq.(42) by MC simulations, as shown in Fig. 3. When At is not too large, 
Eq.@ holds. 

Here we comment on the 6th order integrator. In the MC simulations for Fig. 3, the 
standard construction scheme of Eq.(28) was used. We have other 6th order integrators 



defined by Eq. (|32|) which have less computational cost. These have the same power of 
At to the leading term. The proportional coefficient C n ( here n = 6 ), however, can be 
different for each 6th order integrator. We numerically calculated C n of each integrator 
using quenched QCD and quenched Schwinger ( QED2 ) models. Fig.4 shows (AH 2 ) 1 / 2 
as a function of At, where Stand indicates the standard construction scheme of Eq.(|2£ 



Others are from Eq.(32). The most efficient one is the integrator Yl in Table 1, of which 
C n is roughly a factor of ten smaller than others. The results from others are more or less 
same. The same conclusion is also applied for the Schwinger model as shown in Fig. 5. 
It seems that Yl is always the most efficient one. Thus in the following analysis we use 
Yl as our 6th order integrator. 

Using Eq.(|42|) without higher order terms, we find a formula for the average accep- 
tance, instead of Eq.(|3^), as 

(P acc ) = exp(-C' n uUt"), (43) 
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where C n = C n /v27r. Using this expression, we can easily obtain the optimal step size 
as 

At o P t = r/^-r- (44) 

V nC n V2 

Substituting this result to Eq,(|34"|) we obtain the optimal acceptance as 

(Pacc)o P t = exp(--) (45) 
n 

0.61 2nd 

0.78 4th (46) 
0.85 6th. 



Finally from Eq.(^4|) and Eq.(45) we obtain the optimal efficiency of the n-th order 
integrator: 

( E ff) n opt th = ^M-^{l^-T, (47) 



nC n V2 



and find that Eq.(H^) is governed by only one unknown C n 



The result of Eq.(]45[) suggests that there exists an optimal acceptance depending 
only on the order of the integrator, not on the model. The optimal acceptance increases 
with increase of the order of the integrator. We verify this feature by MC simulations. 
First we show results of QCD for 2nd order one and then that of quenched Schwinger 
model ( QED2 ) for 2nd, 4th and 6th order ones. The Schwinger model is used to reduce 
the CPU time. 

In Fig. 6 results for the 2nd order integrator from three parameter sets among different 
(5 = (5.3, 5.0, 0.0), k = (0, 0.2) and volume V = (6 4 , 4 4 ) are plotted. For all the cases the 
optimal acceptance locates around a value between 60 — 70% which is in good agreement 
with the analytic estimate ( 61% ) of Eq . ([46|) . The results with different lattice volume 
compare lattice size dependence on (Eff) opt . Since QCD is a 4-dimensional model ( 
V = L 4 where L is lattice size ) we find (Eff) op t ~ 1/L for the 2nd order integrator. The 
MC results from 4 and 6 4 lattices with slightly different f3 agree with this expectation, 
i.e. (E ff ) opt [4 4 }/(E ff ) opt [6 4 } ~ 1.5 ( See Fig.6 ). 

Fig. 7 compares (-E77) among 2nd, 4th and 6th integrators. It is clearly seen that 
the optimal acceptance increases with increase of the order, which is again in agreement 
with Eq.( 



5 Comparison of various integrators 

In this section we give a criterion to compare various integrators. Let us consider the n- 
th and m-th order integrators (n > m). As seen in the previous section, each integrator 



has the optimal efficiency given by Eq.(47). In order to have a better performance for the 
n-th order integrator than for the m-th one, the optimal efficiency of the n-th one should 
be larger than that of the m-th one. Furthermore we must consider the cost to perform 
the higher order one since the higher order one needs more arithmetic operations. Thus 
the following equation should be satisfied, 

( E ff)opt th > knm(Eff)™ p J th (48) 



S 



where k nm is a relative cost factor needed to implement the n-th order integrator against 
the m-th one. 

Substituting Eq.([47|) to Eq.(f4|), we obtain 



exp(— ) ? } l > k nm exp(- — ) W — — — r . (49) 
n V nC n V* m V mC m V2 



Rewriting Eq.(|49|), we obtain an expression for volume size with which the n-th order 
integrator performs better than the m-th order one, 

> (k nm eM~- + ~)) nm {-^A (nC n ) m . (50) 
m n \mC m ) 

In the present study we compare (A): 2nd and 4th order integrators, and (B): 4th 
and 6th order integrators. 

(A) : 4th order versus 2nd order 

In this case n = 4 and m = 2. From Eq.(j3l|) we find the relative cost is k±2 = 3. 
Substituting k^p = 3 into Eq.(j5C|) we obtain 

V \ > (3exp(-i)) 4 [l\ C 4 . (51) 

(B) : 6th order versus 4th order 

In this case n = 6 and m = 4. Since we use the scheme of Eq.(^), the relative cost 
/c6,4 is 7/3. Thus we obtain 

Vi >(^M-^Y' '(£)' '(IS,)' '. (52) 
6 Lattice size for higher order integrator 

Now we come to the stage of determination of lattice sizes which are suitable for the higher 
order integrators. Eqs.(|5l|)-(^2|) determine regions where the higher order integrators 
perform better than the lower one. In practice we solve Eqs.(|5T|)-(|52"|) equating both sides 
of the equations. The solutions ( in lattice size ) form a boundary which separates two 
regions: higher order and lower order preferred regions. Unknown C n should be obtained 
from numerical simulations. We choose a small enough At and compute (AH 2 ) 1 / 2 . 
Applying Eq.(||) for (AH 2 ) 1 / 2 we extract C n (= V^rC n ). 

First we study quenched QCD where C n are determined as a function of f3, and then 
go to full QCD with two flavors of Wilson fermions where C n are a function of n. 

6.1 Quenched QCD 

Numerical simulations were performed on a 4 4 lattice. Fig. 8 shows C n as a function of 
(5. We also used an 8 4 lattice at several f3 to check the volume dependence appearing in 
Eq.p^). The values of C n obtained from both lattices were same, which suggests that 
we can get reliable values of C n on the lattice of this size. 
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Fig. 9 shows results of the boundaries determined from Eqs,(51)-(|52[). The circle 
symbols form a boundary which separates the 2nd order preferred region ( lower region 
) and the 4th order preferred one ( upper region ). Similarly the squares separate the 
4th order one ( lower ) and the 6th order one ( upper ). 

The boundary between the 2nd order one and the 4th order one (B2-4) increases as (3 
increases and reaches a plateau at j3 ~ 3.0. The corresponding lattice size is about ~ 10 4 . 
On the other hand the boundary between the 4th and the 6th one (B4-6) decreases with 
(3. At (3 ~ 5.0 the lattice size on B4-6 is about 20 4 . At f3 ~ 5.0, the 4th order integrator 
becomes efficient for L > 10 and the 6th order one for L > 20. 



6.2 Full QCD 

We use a model with two flavors of Wilson fermions. To consider fermion dynamics only 
we take (3 = and simulate the model with varying k. An advantage of taking [3 = is 
that the critical kappa is known analytically, i.e. k c = 0.25. Using the critical kappa the 
quark mass at (3 = is defined by fh q = m q a = ln(l + — 1/« c ))q 

Fig. 10 shows C n as a function of quark mass. The simulations were done on a 4 4 lat- 
tice. We see that C n behaves as a power of quark mass, i.e. C n oc m~ a . Using three data 
at small m q ( K = 0.215, 0.225, 0.230), a are estimated to be a ~ 1.55(8), 4.58(11), 7.65(11) 
for 2nd, 4th and 6th order respectively. For the 2nd order, this result is consistent with 
that of the staggered fermion a 1.5. 

Fig. 11 shows results of B2-4 and B4-6. Both B2-4 and B4-6 increase with de- 
creasing quark mass. We estimate quark mass dependence of the boundaries as: B2-4 
oc C\j 2 ' /C2 ~ m~ a74 and B4-6 oc Cq/C^ 2 ~ m~ ' 73 . For both, values of the power are 
negative, which means that at small quark masses, the lattice size needed to have a gain 
with the higher order integrator increases. If we stay at L < 24(40) which is a lattice size 
available for the current ( or near-future ) full QCD simulations, the 4th order integrator 
can be efficient for fh q > 0.3(0.1). The chiral limit ( m q — > ) is the primary interesting 
case in full QCD simulations. Therefore our results show that the standard leapfrog ( 
2nd order ) integrator is the best one for most full QCD simulations except for heavy 
quarks. 

For comparison we also give results of the Schwinger model with staggered quarks at 
(3 = 0.0. Fig. 12 shows C n as a function of staggered quark mass. The behaviour of C n 
is similar to that of the full QCD case, i.e. C n oc m~ a : a is estimated to be 1.36(9), 
3.36(27) and 5.41(23) for 2nd, 4th and 6th integrators respectively. The estimation is 
based on the data at small quark masses. Fig. 13 shows B2-4 and B4-6. Here note that 
V = L 2 since the Schwinger model is a 2 dimensional model. 

Although we considered (3 = 0.0 case only, at finite (3 we expect the similar diagram 
as in Fig. 11 for small m q a. The reason is the following. The fermionic force is expected 
to become dominant over the gauge force at small m q a. In such a case, the integration 
error from the fermionic part also becomes dominant at small m q a. As seen in Fig. 8 all 
C n of quenched QCD at (3 ~ 5.0 are less than 10. On the other hand all C n at m q a < 0.2 
are bigger than 100. Therefore the naive expectation is that the dominant contribution 
to (AH 2 ) 1 / 2 at small m q a comes from the fermionic part and values of C n at small m q a 

*The alternative definition m q a = — 1/k c ) gives similar results 
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may not drastically change from that of C n at (3 = 0.0. No change of C n results in no 
change of the boundaries ( B2-4 and B4-6 ). 

We might also expect that the multiple time scale method [|j does not work for 
small m q a since it integrates finely only the gauge part and the dominant error from the 
fermionic part remains big at small quark masses. 

To justify the above expectations, we calculate C n at various (3 and illustrate how 
the dominant contribution to (AH 2 ) 1 / 2 comes from the fermionic part. For this purpose 
we choose the Schwinger model with staggered fermions, which requires less CPU time. 
Fig. 14 shows C2 versus m q a at (3 = 0.0, 0.5 and 1.0. At large m q a, the fermionic 
contributions are expected to be small and, in the limit of m q a — > 00 the values of C n 
go to those of quenched ones. In contrast to the situation at large m q a, as m q a goes to 
small quark masses, values of C2 at finite (3 become similar to those of C2 at (3 = 0.0, 
which shows that the fermionic contribution becomes dominant and the values of (3 are 
less important. We observed similar results for the higher order integrators. Therefore 
we deduce the similar result to that at j3 = 0.0. 



7 Summary 

We have investigated higher order integrators for HMC with a criterion which compares 
among various integrators. The criterion is governed by one unknown parameter which 
is easily obtainable from a MC simulation. We have made comparison with quenched 
and full QCD models. 

For quenched QCD the 4th order integrator performs better than the 2nd one for 
a lattice with size L > 10 at (3 ~ 5.0. Of course usually the HMC is not used for 
quenched QCD simulations. If one uses a complicated action which is not implemented 
effectively with a local update algorithm the HMC with a higher order integrator could 
be an efficient algorithm for its simulation. 

For full QCD the higher order integrators can be efficient only for large m q a. For 
instance, on a currently accessible big lattice ( L ~ 24 ) the 4th order one performs better 
than the 2nd order one only for m q a > 0.3, which is out of interest for most full QCD 
simulations. Thus the 2nd order one is the best one for the current full QCD simulations. 

The higher order integrators are turned out to be uninteresting at low (3 ( Fig. 9 ) and 
small m q a ( Fig. 11 ). This may be due to the fact that QCD becomes more perturbative 
(i.e. close to a integrable system) at high (3 and large m q a. The Hamilton's equations 
can be effectively integrated with higher order integrators in perturbative region. On 
the other hand, at low (3 and small m q a (in non-perturbative region), the higher order 
integrators may not be efficient enough. 

The optimal acceptance strongly depends on the order of the integrator. An inter- 
esting case is the 2nd order one, where the optimal acceptance is measured to be around 
60-70% ( analytically estimated to be 61% ). This indicates that a very high acceptance 
like 80-90% or more is not required for the HMC with the 2nd order integrator. We 
suggest to take an acceptance around 60-70% for 2nd order HMC simulations of any 
model. 

We have not considered finite temperature case. Ref.[15] found that quark mass 
dependence of the coefficient ( here ~ C2 ) is very weak in the finite temperature phase. 
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If this is also true for the higher order integrators the boundaries may not increase rapidly 
with decreasing quark mass as fast as in the zero temperature case. 
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APPENDIX 

When the argument x is small, the error function erfc(x) is approximated as follows. 



X 



erfc(x) = 1 — erf(x) 

= 1 — I dte 

\/TT JO 
2 > „2 



1 -=(x + x z + ...) 



7T 

exp( i=x). (53) 

\/7T 



Thus, using (AH) m \(AH 2 ) at small (Ail), Eq.(|37|) is approximated as 

(P acc ) = erfc(i(A J ff) 1 /2) (54) 

« erfc((-A^ 2 ) 1 / 2 ) (55) 
8 

» exp(--l=(A J ff 2 ) 1 /2). (56) 
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Figure 1: Comparison of the average acceptance between exp(— (AH 2 ) 1 / 2 /\/2ti) and 
MC results as a function of (AH 2 ) 1 / 2 . 
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Figure 2: Difference between exp(— (AH 2 ) 1 / 2 / \/2tt) and MC results. Error is denned 
by (MC data - x)/x, where x = ex.p(-(AH 2 ) 1 / 2 /V2tt). 
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Figure 3: (AH 2 ) 1 / 2 for the 2nd, 4th and 6th order integrators as a function of At. The 
simulations ( quenched QCD ) were done on a 4 4 lattice at (5 = 5.0. Three lines in the 
figure are shown to guide the eye. 
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Figure 4: Comparison of (AH 2 ) 1 / 2 among various 6th order integrators as a function 
of At. Stand indicates the standard construction scheme of Eq.(^). Y1-Y3 indicate 
Yoshida's construction scheme. The quenched QCD simulations were done on a 4 4 lattice 
at = 5.0. 
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Figure 5: Same as in Fig.4 but for quenched Schwinger simulations on an 8 2 lattice at 
P= 1.0. 
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Figure 6: The efficiency Eff of the 2nd order integrator as a function of the average 
acceptance. 
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Figure 7: Comparison of the efficiency Eff among the 2nd, 4th and 6th order integrators 
for the quenched Schwinger model as a function of the average acceptance. Simulations 
were done on a 32 2 lattice at (3 = 10.0. 
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Figure 9: Results of the boundaries B2-4 and B4-6 for quenched QCD. Lines are shown 
to guide the eye. 
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Figure 10: C n for full QCD as a function of m q a, where m q a = ln(l + — 1/k c )). 
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Figure 11: Results of the boundaries B2-4 and B4-6 for full QCD at (3 = 0.0. The 
dashed line is an anticipated boundary of B2-4. 
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Figure 12: C n for the Schwinger model with staggered quarks at (3 = 0.0. 
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